function [Axis_rot, Axis_x, Axis_y, Axis_z] = coord_def(point1, point2, point3)
%each point is a three element vector

x = [1 0 0]; y = [0 1 0]; z = [0 0 1];

Axis1 = point1 - point2;
Axis1 = Axis1/norm(Axis1);
Axis2 = point2 - point3;
Axis2 = Axis2./norm(Axis2);
Axis_x = Axis1/norm(Axis1);
Axis_z = cross(Axis1,Axis2);   Axis_z = Axis_z/norm(Axis_z);
Axis_y = cross(Axis_x,Axis_z); Axis_y = Axis_y/norm(Axis_y);

Axis_rot =   [dot(x, Axis_x) dot(x, Axis_y) dot(x, Axis_z) point1(1);...
              dot(y, Axis_x) dot(y, Axis_y) dot(y, Axis_z) point1(2);...
              dot(z, Axis_x) dot(z, Axis_y) dot(z, Axis_z) point1(3);...
              0                 0                   0            1];